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Abstract 

In the presented field line map approach the simulation domain of a tokamak is covered with a cylindrical grid, which is Cartesian 
within poloidal planes. Standard finite-difference methods can be used for the discretisation of perpendicular (w.r.t. magnetic field 
lines) operators. The characteristic flute mode property {k\\ of structures is exploited computationally by a grid sparsification 

in the toroidal direction. A field line following discretisation of parallel operators is then required, which is achieved via a finite 
difference along magnetic field lines. This includes field line tracing and interpolation or integration. The main emphasis of this 
paper is on the discretisation of the parallel diffusion operator. Based on the support operator method a scheme is constructed which 
exhibits only very low numerical perpendicular diffusion. The schemes are implemented in the new code GRILLIX, and extensive 
benchmarks are presented which show the validity of the approach in general and GRILLIX in particular. The main advantage of 
the approach is that it does not rely on field/fiux-aligned, which become singular on the separatrix/X-point. Most tokamaks are 
based on the divertor concept, and the numerical treatment of the separatrix is therefore of importance for simulations of the edge 
and scrape-off layer. 
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1. Introduction 

The modelling of the edge and scrape-off layer of tokamaks 
is in many ways more difficult than the core |[T|. However, this 
region is of high importance since it may have a significant in¬ 
fluence also on the core region, e.g. plasma and impurity den¬ 
sities are often largely set by edge conditions and in important 
operating conditions the edge plays a key role in the improve¬ 
ment of confinement l^. Moreover, a prediction of heat fluxes 
on the divertor plates for future tokamaks is of high importance 
from the engineering point of view l[3l|4]|. 

A major complexity at the modelling is introduced by the 
complex geometry of diverted machines. Field-aligned coordi¬ 
nates are often employed in simulations, since they allow for 
a convenient way to computationally exploit the characteristic 
flute mode property {k\\ ^ kj_) of the structures. However, field- 
aligned coordinates become singular on the separatrix and sim¬ 
ulations cannot span a domain across it. Any set of poloidal (Og) 
and toroidal ((fs) straight field line angles has to satisfy along 
magnetic field lines the condition Q : 


dcps qW' 

At the X-point the poloidal magnetic field vanishes and there¬ 
fore on the separatrix the safety factor q diverges. The straight 
field line angles, which have to satisfy condition (H’ cannot 
span the whole separatrix. As exemplified in fig. ^ the con¬ 
tours of 6s are sucked into the X-point (see also 0). 
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Also often employed are coordinates, where the field- 
alignment property is given up, but which is still aligned 
with flux surfaces, i.e. p(0^) is retained as a radial coordinate. 
However, flux-aligned coordinate systems are still singular on 
points, where Vif/ = 0, i.e. at 0-points and X-points Q. Al¬ 
though these singularities can be cured numerically, O- and 
X-points remain somewhat exceptional points of the numeri¬ 
cal grid (see fig. This could in the worst case even lead to 
numerical artefacts. Moreover, structured flux-aligned meshes 
have a huge resolution imbalance within the poloidal plane due 
to the flux expansion near the X-point El 13. Simulations might 
suffer from this as perpendicular operators arising in practically 
any plasma model (e.g. V^, v^ • V) act approximately isotropi¬ 
cally within poloidal planes of tokamaks. 

The field line map approach is presented in section Al¬ 
though field/fiux-aligned coordinates may become singular, the 
operators appearing in plasma models are still well defined, of 
course. The idea behind the approach is that the flute mode 
property can also be exploited at the discretisation step without 
any need for construction of a field/fiux-aligned coordinate sys¬ 
tem. The approach consists of a cylindrical or Cartesian grid 
with a field line following discretisation for parallel operators. 
A separatrix can be treated as well as a magnetic axis, where 
X/O-points are treated like any other grid point and no reso¬ 
lution imbalance arises. Ultimately, the result is similar to the 
flux-coordinate independent (FCI) l[T0l[TT][T2l approach. How¬ 
ever, as the motivation for the FCI approach was initially still 
based on field-aligned coordinate systems, the derivation is here 
performed completely without any reference to field- or flux- 
aligned coordinate systems. 

As the discretisation of perpendicular operators in the field 
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Figure 1: a) Contours of flux label p (black) and poloidal straight fleld line angle 
9s (red, symmetry coordinates) in diverted geometry for equilibrium from QS). 
b) Example for flux aligned mesh. The contours of the poloidal coordinate are 
orthogonal to flux surfaces. The X-point is connected to eight cells instead of 
the usual four cells. 


line map approach is straight forward, the main emphasis in this 
paper is on the discretisation of parallel operators, especially the 
parallel diffusion operator in section [2^ Since an interpolation 
or integration is involved at the discretisation, parallel operators 
exhibit also numerical perpendicular ’diffusion’. Motivated by 
previous work from ifT^fT^ . a numerical scheme is developed 
which exhibits very low numerical diffusion. The discussion 
extends previous work from ca. Several model problems are 
also discussed in [Appendix A| 

The developed numerical methods are implemented in the 
new code GRILLIX. In section extensive benchmarks per¬ 
formed with GRILLIX are presented, which show the validity 
of the field line map approach in general and GRILLIX in par¬ 
ticular. 

The paper is concluded with a summary and final remarks in 
section [H 

2. Field line map approach 

2.7. Overview 

The field line map approach is described in the following for 
the case of a toroidal configuration (/?, Z, ^), but it can be ap¬ 
plied also to axial periodic configurations (v,y,z), where z is 
the axial coordinate. The transition should be trivial. 

For a tokamak a cylindrical coordinate system is well de¬ 
fined everywhere, except for the toroidal symmetry axis which 
is outside the domain of interest. We span the simulation do¬ 
main with a cylindrical grid Rp Zj, (f^ (For axial configurations 
a Cartesian grid Xi,yj,Zk is used). Within each poloidal plane 
k the grid [Ri,Z^ is Cartesian and bounded by extreme flux 
surfaces, which is the only dependence on flux surfaces of the 
approach. Based on the assumption of a strong toroidal field 
^^tor ^ any perpendicular operator is approximated by 

derivatives with respect to only R and Z, but not ip. In order 


to exploit the flute mode property k\\ a dense resolution 

is chosen within poloidal planes, whereas the grid is sparsified 
along the ip direction. The low resolution in ip requires a field 
line following discretisation for parallel operators to achieve a 
sufficient directional accuracy for parallel operators. 

It is noted that the assumption of a strong toroidal field is 
not a strictly necessary condition. If this assumption breaks 
down, a dense resolution also in ip would have to be retained, 
and the perpendicular operators would have to be adjusted to 
take into account also derivatives with respect to ip. Ultimately, 
the field line map approach would go over into a discretisa¬ 
tion on a dense cylindrical grid, where the flute mode property 
would not be exploited any more. However, the method would 
still retain its validity. 

2.2. Perpendicular operators 

Under the assumption of a strong toroidal field, any perpen¬ 
dicular operator can be approximated with derivatives with re¬ 
spect to only R and Z, e.g. the perpendicular Laplace operator 
becomes: 
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where the inverse metric of the cylindrical coordinate sys¬ 
tem and If the contravariant components of the unit vector of 
the magnetic field. Hence, the Stencil of any perpendicular op¬ 
erator remains within the Cartesian poloidal planes. Standard 
finite difference methods can be used for their discretisation 
(e.g. ifTTll ). e.g. a second order finite difference method yields 
for the discrete perpendicular Laplace operator : 
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where h denotes the Cartesian grid spacing. Bold face is 
used on the discrete level for vectors and matrices, i.e. u := 
(i^i,i,i, ^2,1,1, • • • the subscript i,j,k denotes the corre¬ 

sponding grid point. 


2.3. Parallel gradient 
The parallel gradient is: 


V\\u = lim 
" 6^0 


u (x^) - u (x) 

6 


( 4 ) 


where e is the arc length along a magnetic field line and: 
dx^ 

=b, x^(0) = X. (5) 

de 

This motivates the discretisation of the parallel gradient via a fi¬ 
nite difference along magnetic field lines. The Stencil will cover 
neighbouring poloidal planes, and since the grid is sparsified 
along the toroidal direction, a field line following discretisation 
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has to be performed. Field lines are traced from each grid point 
towards the neighbouring poloidal planes according to: 

+A^ +A^ 

0 0 

with the contravariant components of the magnetic 

field. The integrations start from grid points and '+' 

denotes co/counter-direction with respect to magnetic field. For 
axisymmetric equilibria the penetration points Zfj) are in¬ 
dependent of the toroidal grid index k. Also the lengths along 
field lines are computed according to: 
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where B^^^ = ^JWB^ = RB"^ the toroidal field strength. 

A sketch of the discretisation of the parallel gradient is shown 
in fig. 1^. The discrete parallel gradient is collocated at posi¬ 
tions half way along the magnetic field line towards the neigh¬ 
bouring poloidal planes. This gives rise to two possible dis¬ 
cretisations {'V and Since the penetration points do not 
in general coincide with grid points, a 2D-interpolation has to 
be performed within the Cartesian grid, to obtain the values for 
some quantity u at the penetration points. The interpolated val¬ 
ues are denoted with The discrete parallel gradient opera¬ 
tors Q- are defined as: 


, ik ^Uj,k 

(Q-u). .. := +— - - 




( 8 ) 


The parallel gradient at the considered grid point itself could be 
obtained via a further linear interpolation of the discrete paral¬ 
lel gradient along the magnetic field line. However, since this 
work will concentrate on the discretisation of the parallel diffu¬ 
sion operator, this issue is left for future work. Results on the 
parallel gradient can also be found in lfTT]| . 

Later on, especially for field lines which are strongly dis¬ 
torted, another discretisation of the parallel gradient is useful, 
which is based on integration: 

1 1 1 r 

V\\u = —V • (uB) = — lirn — uB - dS. (9) 

dv 


The surface integration can be mimicked on the discrete level. 
Flux boxes around magnetic field lines are taken as discrete fi¬ 
nite volumes (see fig.[^), such that the only contributions to the 
surface integral come from the toroidal ends of the flux box: 
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where are the discrete volumes, is the magnetic field 
strength in the center of the flux box and AA -. is the surface 
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Figure 2: Sketch of discrete parallel gradient: a) Interpolation method, b) inte¬ 
gration method. 


overlap of grid point (n,m,k± 1) with the toroidal end of the 
flux box surface of grid point (/, j, k) (Toroidal index k can be 
dropped in due to axisymmetry). Based on the fact 

V • B = 0, the flux box volumes can be computed with 
high accuracy during the field line tracing process. 

The introduction of discrete fluxes i p ^2 11 • • •) 

useful for later purposes: 


Q-u, 


( 11 ) 


where Q- are matrices according to equation ([^ or (10), which 
contain also the interpolating coefficients or the surface over¬ 
laps. 


2.4. Parallel diffusion 

In this section the discretisation of the parallel diffusion op¬ 
erator, defined as 


D\\u \= V • [bV||w], (12) 


is discussed. 

The interpolation or integration process at the discretisation 
of the parallel gradient introduces an erroneous numerical per¬ 
pendicular coupling among distinct field lines. Therefore, the 
discrete parallel diffusion operator will exhibit a spurious nu¬ 
merical perpendicular ’diffusion’, which depends in general on 
resolution and the applied interpolation/integration method. In 
a magnetised plasma the dynamics can be strongly anisotropic. 
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i.e. very fast along magnetic field lines and very slow perpen¬ 
dicular. The numerical perpendicular diffusion arising from the 
discretisation error of the parallel diffusion operator might then 
compete or even overwhelm the real perpendicular diffusion 
or dynamics. Therefore, it is important to construct numeri¬ 
cal schemes which exhibit only very low numerical perpendic¬ 
ular diffusion. In the following, two schemes for the discrete 
parallel diffusion operator are presented, a naive discretisation 
and a discretisation according to the support operator method 
|[T8l[T3|. The support scheme was motivated from (131 [Ml, and 
it will turn out that it exhibits a much lower numerical diffu¬ 
sion. In many cases, where an increase of resolution might not 
be possible due to hardware constraints, the naive discretisation 
might fail, whereas the support schemes still work well. 


2.4.1. Naive discretisation 

Under the assumption V • b ky, the parallel diffusion oper¬ 
ator can be approximated as: 

£>II - V2. (13) 


This motivates a discretisation via a further finite difference of 
the parallel gradient along magnetic field lines. The discrete 
parallel diffusion operator is proposed as: 
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2.4.2. Support operator method 

The naive scheme yields a consistent discretisation, but it 
does not conserve any ’good’ property of the parallel diffu¬ 
sion operator. Based on the support operator method |[T^IT9ll 
a scheme is derived which conserves the self-adjointness of the 
parallel diffusion operator on the discrete level. 

Let u, V be two real valued arbitrary scalar fields, we define a 
scalar product via an integration over the whole domain: 

<M,v> := / uvdV. (15) 

V 


It holds that: 


{u,D\\v) = J" t/V • [V||vb] dV = J" uVjjvh-dS-J" V||t/V||vJy 

V 

(16) 


dV 


If we further assume that the domain is periodic (^-direction) 
and/or the quantities vanish at the boundaries, the following in¬ 
tegral equality holds: 

<t/,D||v) = -<V||W, V||v), (17) 


+ + + + 



Figure 3: Illustration of discrete volumes. AV^ are flux boxes 
toroidally limited by \^ipk{A) - A(pl2,(p]^(^x} + and toroidally limited 

±Aip\ 


is discretised, and derived operators are constructed via a dis¬ 
crete analogue of certain integral equalities. In our case we will 
derive the discrete parallel diffusion operator with the parallel 
gradient as the prime operator and via conservation of integral 
equality 'GZl' on the discrete level. We discuss the construction 
of the discrete parallel diffusion operator for the case of of the 
'-f' discretisation in the following. The discretisation fol¬ 
lows analogously and a combination of both is given at the end 
of this section. 

Let us start from the discrete parallel gradient according to 
equation ^ or ( [TQ| ), which can both be represented with a ma¬ 
trix . Scalar functions u, v are collocated on the basic cylin¬ 
drical grid. Gradients q^, are collocated on points half way 
along magnetic field lines towards neighbouring poloidal planes 
(see again fig.[^. The discrete parallel gradient maps from the 
cylindrical grid (SG) to that gradient’s grid (FG^): 

Q^:SG^FG^ (19) 

In the next step we want to mimic the integral equality ( [TT] ) on 
the discrete level. However, on the discrete level the left hand 
side of equation ( [TT] ) denotes a scalar product on the space S G, 
whereas the right hand side on the space FG^. Therefore, we 
have to define two discrete scalar products, which both mimic 
an integration over the same whole domain: 

vsG ■= L “-I w := L ’ 

A fi 

( 20 ) 

where Greek letters indicate summations over all grid points. 
The finite volumes are chosen as finite flux boxes around mag¬ 
netic field lines as illustrated in fig. The discrete parallel 
diffusion operator. 


It is immediately obvious that: 

Vj = -V.[bo], Dj=D||. (18) 

The method of support operators gives a procedure for the con¬ 
struction of discrete second order operators which conserve cer¬ 
tain integral equalities on the discrete level. Within the sup¬ 
port operator method a first order operator, the prime operator. 


^supp,+ . sg^SG, (21) 

can now be derived from the discrete parallel gradient, and by 
imposing the integral equality 'GZl' on the discrete level. On the 
one hand: 

A,cr 
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On the other hand: 


- <Q+u, (23) 

H,V,T 


The conservation of integral equality (p^ on the discrete level 
demands the equality of expressions (| 22 |) and for arbitrary 
u, V. Via relabelling of the indices v ^ A, r ^ cr the discrete 
parallel diffusion operator can finally be obtained: 


YySUpp,+ 
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In axial geometry the volume factors vanish 
(aVi = A^^ = h^Az^ and the discrete parallel diffusion 
operator becomes: 


^ _^q+^7-q+^ (25) 


two-dimensional model problem is considered. The coordinate 
V plays the role of a coordinate within poloidal planes and z the 
role of a periodic toroidal/axial coordinate, and the domain is 
spanned by a regular grid v/, Zk, with grid spacing h in the v- 
direction and grid spacing Az in the z-direction. The magnetic 
field is uniform with an inclination with respect to the grid: 

b = sin 6ex + cos (29) 

The displacement of the penetration points with respect to the 
grid point can be expressed by a factor /: 

xf = Xi ± fh, / = ^ tan 6, (30) 

Without loss of generality only cases of slight inclination with 
respect to the grid are considered in the following, i.e. |/| < 1 . 
Using a linear interpolation the values at the penetration points 
are determined as: 


which expresses the self-adjointness property on the discrete 
level. For toroidal geometries, due to the general geometry and 
the l/R dependence of the toroidal field strength, the volume 
factors account for effects arising from V • b 9 ^ 0. In addition 
the volume factors matter if the parallel distances Asfj vary, 
i.e. if the grid is not equidistant in the parallel sense. 

Finally, a small modification is applied to the scheme. A dis¬ 
crete parallel diffusion operator can also be derived analogously 
with the discretisation of the parallel gradient. This will end 
up in a generally different but consistent scheme by itself. How¬ 
ever, it is desirable that the final scheme is independent of this 
initial arbitrary choice. To achieve this, the average between 
both schemes is taken: 

;= 1 . (26) 

This modification does not alter the self-adjointness property 
on the discrete level, since: 


(u, D7^v\ = - (u, 07 '’’^ v) + - (u, D7^’-v\ 

\ ’ II /SG 2 \ II 'SG 2 \ II fSG 

= f(D7'’’Vv\ (d7'’’'u,v\ =(d7Vv) . 

9 \ II /SG 9 \ II /SG \ II ’ /SG 


SG 

(27) 


A further remark concerns numerical stability. The support 
operator discretisation excludes possible numerical instabilities 
arising from the interpolation, since it guarantees a strict de¬ 
crease of the L? norm: 


Kk = (1 “ f)^i,k±l -^f^i±l,k±G (31) 


The parallel distances are all equal A^-^ = A^- = Az^, 

as also the finite volumes AVa. = A'V^ = hAz. To illustrate 
also the construction of the scheme we give for this case also 
explicit expressions for the matrices Q- on a 3 x 3 grid: 

/ -1 1-/ / \ 

-1 1-/ / 

-1 1-/ 

-1 1-/ / 

-1 1-/ / , 

-1 i-f 

1-/ / -1 

1-/ / -1 

1-/ -1 / 

(32) 

The only inner grid point is / = 2, k = 2, which corresponds to 
the fifth row of the matrices Q-. The discrete parallel diffusion 
operator at the inner grid point is illustrated in fig. 1 ^ for the 
naive scheme according to equation ( p^ a nd in fig. for the 
support scheme according to equation (|26|). It is apparent that 
the Stencil of the support scheme is larger. In the limit that 
the penetration points coincide with grid points, i.e. / = 0 , 1 , 
both schemes yield the standard second order finite difference 
expression. 

Now we investigate the action of the discrete parallel diffu¬ 
sion operators on a Fourier mode u = exp (ikxX ik^z). The 
analytic result is: 

D\\u = with: k\\ := kx sin^ -r cos 6. (33) 




f <Q"u,Q u)to- < 0, 
(28) 

which is in generally not fulfilled with the naive discretisation 
method. 


Performing a Taylor expansion in (kxK k^Az) yields: 
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U, (34) 
(35) 


2.5. A two-dimensional model problem 

To illustrate the difference between the naive and the support 
scheme, and to enlighten the advantage of the support scheme a 


It is apparent that for finite displacement (/ 9 ^ 0,1) the lead¬ 
ing error for the naive scheme scales like {kxhf /As^. This re¬ 
flects the error of the linear interpolation, which scales also like 
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Figure 4: Model problem for homogeneous magnetic field. Discrete parallel 
diffusion operator for a) naive scheme, b) for support scheme. 


on, we consider some blob with finite poloidal extent of order 
h at some plane. In reality the blob should diffuse along mag¬ 
netic field lines and thus spread across many grid points of the 
neighbouring poloidal plane. For the naive scheme this does 
not cause any problems, since the discrete parallel diffusion 
operator is computed by just ’taking’ values from neighbour¬ 
ing poloidal planes. However, in the support scheme a value is 
not only ’taken’ but also ’sent’ towards neighbouring poloidal 
planes. In other words, if the parallel gradient is computed via 
interpolation at the penetration points, grid points which lie in 
between might not be connected to the original points by the 
scheme. The blob does not spread properly across the grid 
points, but only diffuses to points which are connected by the 
scheme. Finally, erroneous wiggles arise in the neighbouring 
poloidal plane. 




(kxh)^. The error term represents a numerical perpendicular dif¬ 
fusion with diffusivity coefficient: 

naive _ /(^ ~ /)^ 

A ±,num ^ ^ 


Figure 5: a) Effects of map distortion on numerical scheme. In the illustration 
only thick points are connected by the support scheme (linear interpolation). 
A blob should in reality spread smoothly across many grid points. However 
the support scheme produces unphysical wiggles (dashed) in the neighbouring 
poloidal plane, b) Examples for conformal distortion (top) and angular distor¬ 
tion (bottom). Grey squares indicate base cells within the neighbouring poloidal 
plane. 


For the support scheme an improved scaling of 
(kxh, kzAz)^/As^ for the leading error is obtained. Therefore, 
the numerical diffusion becomes actually a hyperdiffusion. 

The same analysis performed here for linear interpolation is 
repeated in Appendix A.3 with third order polynomial interpo¬ 
lation. Moreover, in Appendix A.l a model problem for the 
case of an inhomogeneous magnetic field is discussed, and in 
[Appendix A.2| an example for the application of the integration 
method is given. 


2.6. Map distortion 

Especially close to the X-point field lines become strongly 
distorted. To illustrate this, we consider the lowest order ex¬ 
pansion for the magnetic field around the X-point 1201 : 

B = Bo [cj + O' {xe^ - ye,,)], (37) 

Two generic field lines, which have initially a poloidal distance 
do, will exponentially diverge: 


We define a criterion in order to exclude numerical artefacts 
arising from the map distortion: Starting with an initial square 
of lateral length h at some plane, the square encounters a dis¬ 
tortion, as its edges are traced. Therefore, also the map of 
the square at the neighbouring poloidal plane is distorted (see 
fig .|^). Neglecting the weak variation of the area of the 
square is thereby conserved due to flux conservation. We can 
then distinguish two types of distortion 1 ^ . a conformal where 
the square is stretched in one direction and squeezed in the other 
resulting in a quad with disparate lateral lengths, and an angular 
distortion where e.g. two angles become acute and the other two 
obtuse resulting in a parallelogram. We quantify the distortion 
by: 


longest side of mapped quad /, j 
ij shortest side of mapped quad /, j ’ 
largest angle of mapped quad /, j 
ij smallest angle of mapped quad /, j 


(39) 

(40) 


d(z) ~ doexp(az) (38) 

We clarify now the effect of such a strong distortion on the nu¬ 
merical scheme. 

We consider again a two-dimensional setup as illustrated in 
The penetration points of two neighbouring field lines 
separate exponentially according to expression ( [38] ). Further 


The distortion decreases as the toroidal/axial resolution is 
increased, and a reasonable criterion is to require enough 
toroidal/axial resolution such that the mapped quad does not 
cover more than two squares in each direction within the neigh¬ 
bouring poloidal plane. Otherwise the map might jump across 
grid points. Applying this constraint, results in a threshold for 
the distortion (see fig.[^) of ^ 4 and da <3, but we require 
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for convenience the same threshold for both, i.e.: 


dc, da < 4. 


(41) 


Another numerical artefact may arise, if for two neighbour¬ 
ing (perpendicular) points there is a jump in the points which 


are involved in the interpolation (see example in Appendix A. 1 
around grid point i = I and equation ( |A.7| )). This may cause 
a small oscillation on the grid scale, which might be cured by 
adding a small amount of high order perpendicular diffusion. 

Any problems related with a distortion of the map can be 
cured by using the integration method for the discrete parallel 
gradient (equation (p^). With this scheme the information is 
spread properly across grid points in the neighbouring planes. 
In the illustration of fig.[^ the mapped area is obtained by trac¬ 
ing only the four corners of the base square resulting in a quad 
as mapped area. However, if the map distortion becomes really 
strong even this might not suffice and a more detailed contour 
of the base square might have to be traced which would yield 
a polygon as mapped area. Problems related with a jump of 
the points which are involved in the interpolation are also re¬ 
solved with the integration method (see [Appendix A.2| ). From 
a practical point of view the integration method might seem to 
be second choice compared to interpolation methods, since its 
implementation appears cumbersome though possible as will 
be shown with an example in section |T2T] (figure [T^). Never¬ 
theless, it is presented here to show that even complicated cases 
with strong map distortion do not pose problems of principle to 
the support scheme. Moreover, there might be cases where its 
application might be necessary, if no sufficient toroidal resolu¬ 
tion can be supplied to bring the distortion below the requested 
threshold given in ( |4T] ). 


3. Benchmarks and examples 

The numerical methods presented in section are imple¬ 
mented in the new code GRILLIX. In this section we present 
benchmarks, which shall show the validity of the field line map 
approach in general and GRILLIX in particular. As a model 
problem the parallel diffusion equation 


• S-C: Support scheme with integration method for parallel 
gradient 

Each bilinear interpolation involves four grid points and each 
third order bipolynomial interpolation 16 grid points symmet¬ 
rically arranged around the considered penetration point (see 
e.g. 1221 ). For the S-C method the overlaps AAfj^^ arising in 
the definition of the discrete parallel gradient (expression ( p^ ) 
are computed via the routine given in 1231 . 

3.1. Axial circular equilibria 

For axial circular equilibria the problem ( |4^ can be solved 
analytically. We will consider in the following flux shells, 
i.e. radially bounded domains p g [pmm^Pmax], with p \= 
a flux label. Without loss of generality, we consider 
only a single radial, poloidal and axial mode (r, m, n) as initial 
state: 


u{t = 0, x) = sin nr 


P Pmin \ . 

_ I Cl 


Pmax Pmin 


sin {mO + nz ), 


(43) 


with tmiO = ylx being the poloidal angle. The analytic solution 
is: 


u{t, x) = u{t = 0, x) exp 


All 


m -\- n q(p) 

^(p)2 -rp2’ 


t 

im,n(P^ 


(44) 


with q(p) = B^/B^ the safety factor. A characteristic time 
for the decay of the mode is given by tm,n(po) with po := 
(Pmax Pmin^ 1'^' 

The benchmarks, presented in the following, were performed 
with pmin = 0.1 and pmax = 0.2. Since the goal is the investiga¬ 
tion of the spatial discretisation error, time steps were chosen 
sufficiently small, such that the temporal discretisation error 
played always only a subdominant role. The numerical error 
will be quantified in the norm: 


eiit) \= 


analytic ( 0 - ^numeric mi 
\\^analyticif )\\2 


(45) 


ju=x\\^\u (42) 

is considered. Space scales are normalised to Rq in toroidal ge¬ 
ometries respectively Laxl^n in axial geometries, where Lax is 
the axial periodicity length. Time is measured in R^lxw respec- 
tively LlJ {47:X)- 

Five possible discretisation schemes for the parallel diffusion 
operator are investigated: 

• N-1: Naive scheme with bilinear interpolation 

• N-3: Naive scheme with third order bipolynomial interpo¬ 
lation 

• S-3: Support scheme with bilinear interpolation 

• S-3: Support scheme with third order bipolynomial inter¬ 
polation 


where the norm is computed on the discrete level as \\u \\2 := 
^J{u, u)sG according to equation (20). 


3.1.1. k\\ 9 ^ 0 modes 

Firstly, we consider modes which have a non-vanishing par¬ 
allel wavevector k\\ oc m + n q{p). 

In the case of a pure axial magnetic field {q oo) the pen¬ 
etration points coincide with grid points and the interpolation 
becomes exact. The discrete parallel diffusion operator reduces 
to the standard second order finite-difference expression for all 
schemes. Hence, this benchmark serves merely as a first basic 
test on the correctness of the implementation of the schemes in 
GRILLIX. Fig. shows that the expected second order accu¬ 
racy with axial (^parallel) resolution was obtained with GRIL¬ 
LIX. 

If a finite safety factor q oo considered, the penetration 
points do in general not coincide with grid points any more, and 
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Figure 6: Numerical error (for all schemes exactly overlapped) in dependence 
on axial resolution for a(r= l,m = 0, n = 1) mode with q = oo. 


the numerical error is dependent on the poloidal resolution h 
and the scheme. In fig. the difference between the numeric so¬ 
lution and the corresponding analytic solution is shown for var¬ 
ious different schemes. The basic behaviour is consistent with 


the predictions from the two-dimensional model of section 2.5 


and [Appendix A.3[ The overall error level is largest for the N- 
1 scheme of order O (see equation (34)), whereas the error 

level is lower of order O (h^) for N-3 (see equation (A. 17)), S-1 
(see equation ( [^ ) and S-3 (see equation ( A.18| )). ^Significant 
errors arise at the radial boundaries of the domain. Though the 
parallel diffusion equation ( |42| ) does not require radial boundary 
conditions, some assumption about the quantity must be made 
also outside the domain for the interpolation. Within GRILLIX 
it is implicitly assumed that u = 0 outside the domain. This re¬ 
sults in radially discontinuous derivatives at the boundary, and 
therefore especially for higher order interpolations large errors 
arise close to it. 



Figure 7: Difference between numeric and analytic solution for a 
(r = I, m = 3, n = 1) mode with ^ = 3.4 at r = fsj (for definition of rsj 
see equation j44) ) on plane z = 0 for various schemes. Resolution was 
h = 6-l0-\ Az/27i= 1/32. 


The numerical error in dependence on axial resolution for 


two fixed poloidal resolutions is shown in fig. For low ax¬ 
ial resolutions the error follows a Az^ law showing the sec¬ 
ond order accuracy of the schemes along magnetic field lines 
(Az ~ A^-). For higher axial resolutions the error is dominated 
by the poloidal resolution h and deviates from the Az^ line but 
increases. Again in agreement with the scaling derived from 
the two-dimensional model from section |2.5| and jAppendi}^ 
|A.3| the deviation occurs first for the N-1 scheme and occurs 
later for the other schemes. In fig. [^the error in dependence 
on the poloidal resolution for two fixed axial resolutions is 
shown. Again the convergence is slowest for the N-1 scheme, 
whereas the others exhibit similar convergence rates. Note that 
though the S-1 scheme is based only on bilinear interpolation, 
it seems to perform even slightly better than the N-3 scheme. 
All schemes converge to the same value of ^2 ~ 1-1 • 10“^ at 
Az/27t = 1/32 (fig.[^) and ^2 ~ 3.5 • at Az/Itt = 1/64 
(fig#). These values are settled by axial resolution (see red 
stars in fig.[^). 


a) 



b) 



Figure 8: Numerical error in dependence on axial resolution for a 
(r = 1, m = 3, n = 1) mode with q = 3.4 for fixed poloidal resolution of a) 
h = ^ ■ 10“^ and b) h = 4 ■ 10“^. Red stars indicate values obtained in the limit 
h ^ 0 (see fig.j^. 


3.1.2. k\\ = 0 structures 

We consider an initial state which is Gaussian in the poloidal 
plane and a delta function in the axial direction, i.e.: 


u{t = 0, x) = exp - (x - jwl + {y- ygf jw: 


Hz), (46) 


where the delta function is modelled on the discrete level with 
a Kronecker delta 6{z) Sok^njAz. This structure diffuses 
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a) 


a) 



b) 



Figure 9: Numerical error in dependence on poloidal resolution for a 
(r = I, m = 3, n = 1) mode with q = 3A for fixed axial resolution of a) 
Az/27r = 1/32 and b) Az/2n = 1/64. 


along magnetic field lines until it has established a state with 
^11 = 0, which should remain stable. In fig.H^the temporal evo¬ 


lution of such a Gaussian blob is shown computed with the N-1 
and S-1 scheme. The example illustrates that after the structure 
has established a state in which k\\ = 0 ait ^ 30, it remains sta¬ 
ble for a long time with the support operator scheme, whereas 
it decays fast with the naive scheme due to numerical perpen¬ 
dicular diffusion. 

Any decay of zonal modes (r, m = 0, n = 0) arises from er¬ 
roneous numerical perpendicular diffusion. The numerical dif¬ 
fusion can be quantified by estimating the decay rate of the L^- 
norm of zonal modes. In fig.pT^ the measured decay rate in de¬ 
pendence on the poloidal resolution of the zonal modes is plot¬ 
ted for fixed axial resolution. In fig. \TT\ ) the numerical decay 
rate in dependence on the axial resolution for a fixed poloidal 
resolution of a zonal mode is plotted. Summarizing the numer¬ 
ical decay exponent scales like: 


Ynum ' 


(kphf /Az^, 

[kphf 

kp^ I 


for N-1, 
for N-3, 


for S-1, 

A (kp}^ + B (kph^ 


(47) 


/Az^, for S-3. 


These results confirm again the scalings of the two- 
dimensional model (see equations ( [34| ), ( [35] ), ( A.17 ), ( |A.18 )). 


Only for the S-3 scheme a break in the scaling is observed, 
which arises from the boundaries, where the conditions for 
the interpolation become worse (It has also been observed that 
the numerical decay sets first in near the boundaries). The 




- 6.2 - 0.1 0 0.1 0.2 ■°- B .2 - 0.1 0 0.1 


-0.3 -0.2 -0.1 0 0.1 0.2 0.3 


b) 




S- 1 ,t =30 

• 



• 



-0.3 -0.2 -0.1 


0.1 0.2 0.3 


Figure 10: Snapshots of diffusion of blob at f = 30 (left) and t = 200 (right) 
on plane z = 0 simulated with a) N-1 scheme and b) S-1 scheme. A radially 
constant safety factor of q = 3 was chosen. Resolution was h = 2 ■ 10“^, 
Az/2n = 1/8. Initial blob was located at the bottom (xg = 0, = -0.15, = 

Wy = 0.025 according to equation (j^). 


performed numerical measurements and the discussion of the 


two-dimensional model problem (see section [23] and |Appendix 
|A.3| ) may suggest the general rule that, if p is the order of the 
polynomial interpolation, the numerical decay exponent scales 
like: 


(kphy /Az^, for naive schemes. 


ynum ^ 


(kph^ ^ /Az^, for support schemes. 


(48) 


apart from effects introduced by the boundary, which introduce 
a break into the scaling at high poloidal resolutions. The inves¬ 
tigation of numerical perpendicular diffusion for other interpo¬ 
lation techniques is left for future work. 

3.2. Realistic toroidal geometry 

We consider an equilibrium, which is an analytic solution to 
the Grad-Shafranov equation ifTbl . In order to prove the appli¬ 
cability of the scheme also to complicated and realistic tokamak 
geometries, we consider a flux shell in the edge {pmin = 0 - 90 , 
= 0 . 95 , the normalized radial flux label is defined as 
p \= xjii!/ - iks) / (<Ao - where ij/Q^ the poloidal magnetic 
flux at the magnetic axis, respectively at the separatrix). 
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a) 


a) 


b) 




Figure 11: Numerical decay rate in dependence on a) poloidal resolution of 
mode for fixed axial resolution Az = 27r/2. kp := nr Up max - Pmin) is the radial 
wavevector. For the establishment of this plot the radial mode number r has 
been varied at two distinct poloidal resolutions ofh = 6- 10“^ and h = 4 - 10“^ 
have been used, b) Numerical decay rate in dependence on axial resolution for 
fixed poloidal resolution of mode (kph = 0.75 


3.2.1. Verification of map distortion 

To illustrate the effects of the map distortion on the schemes, 
we consider the diffusion of a Gaussian blob situated at the bot¬ 
tom of the flux shell in proximity to the X-point. In fig. [T^,b 
snapshots are shown of solutions computed with the naive 
scheme N-3 and the support scheme S-3. The simulations were 
performed only with two poloidal planes A^/27r = 1/2. On 
a coarse scale the blob appears at similar positions, but on a 
finer scale strong unphysical wiggles arise in the simulation per¬ 
formed with the support scheme. As discussed in section |2.6[ 
the erroneous wiggles are a result of the strong map distortion, 
which was dp = 44 and da = 138. These problems can be cured 
by requiring enough toroidal resolution, as shown in fig. [T^ 
where the simulation was repeated with a higher toroidal reso¬ 
lution of bxiplln - 20. The distortion remained with values of 
dp - 1.7 and da - 1.9 well below the proposed threshold of 
4 (see expressions (41)). The problems with the map distortion 
can also be cured via the integration method for the parallel gra¬ 
dient (equation (1^)). In fig. 12 i the same simulation was per¬ 
formed with again a low toroidal resolution of bxiptln - 1/2, 
but with the S-C scheme. No erroneous wiggles but smooth 
structures arise. 


The strongest map distortion occurs around the X-point. For 
the given equilibrium the mapped quads around the X-point for 
different toroidal resolutions are illustrated in fig.[T^. For too 



c) d) 



R R 


-0.5 -0.4 -0.3 


0 


0.2 0.3 0.4 0.5 


Figure 12: Snapshots of diffusion of Gaussian blob at r = 10, = tt. Initial state 

was Gaussian blob = 0.88, Zg - -0.29, wr-wz-^- 10“^) and poloidal 
resolution was h = I • 10“^. a) Computed with N-3 scheme and tSpijln - 
1/2, b) with S-3 scheme and tS^lln - 1/2. The strong map distortion causes 
unphysical wiggles. This can be cured by increasing toroidal resolution as was 
done in c) (S-3 scheme with ^Lpjln = 1/20) or d) by using the integration 
method for the parallel gradient (S-C scheme with iSipIln = 1/2). 


low toroidal resolutions (blue) the mapped quads are strongly 
distorted and extend over many grid points. The quality of the 
grid would be much too poor for support schemes (interpola¬ 
tion based, i.e. S-1, S-3). At higher toroidal resolutions (green) 
the mapped quads are only weakly distorted and the support 
schemes would work well. In fig. [T^ the map distortion in 
dependence on toroidal resolution is plotted. A slightly higher 
toroidal resolution than A^/27r =1/16 would have to be used, 
to fulfil the proposed requirement 


3.2.2. Convergence behaviour 

For realistic geometry rigorous convergence checks -as it was 
done for axial circular geometry in section [3TT| are elaborate 
and computationally quite costly. Effects of a complex mag¬ 
netic field structure on the schemes (map distortion) have been 
discussed in detail in section [Z61 and have been verified in sec¬ 
tion 3.2.1 Therefore, only a rough convergence test for the 
schemes is presented here for realistic geometry. 

For the diffusion of the Gaussian blob, which was presented 


10 





















































































































a) 



b) 



factory with the available resolution, whereas N-3, S-1 and S-3 
exhibit similar convergence rates. The naive scheme conver¬ 
gences against a slightly different value 0.42) than the sup¬ 
port schemes 0.40). The reason for this is that also slightly 
different operators are discretised with the naive scheme 
and the support schemes (V • [bV||o]). 


a) 




Figure 13: a) Mapped quads (’+’-direction) for some sample grid points around 
X-point (red cross) for different toroidal resolutions, b) Map distortion around 
X-point in dependence on toroidal resolution. 


in section [TTT| a spatial resolution scan was performed. How¬ 
ever, we only consider the obtained maximum value on plane 
(f = 71 at time t = 10.0, i.e.: \\u(t = 10,R,Z,(p = 7r)||oo. Effects 
of numerical diffusion will manifest themselves in a drop of this 
value. A reduction of the timestep below values where the tem¬ 
poral discretisation error would be negligible turned also out to 
be computationally to costly. Therefore, the timestep was held 
constant at At = 1 • 10“^ regardless of spatial resolution, in order 
to obtain the behaviour of the spatial discretisation error. 


In fig. [T^ the convergence behaviour in dependence on 
toroidal resolution is shown. At low toroidal resolutions the 
support and the naive schemes yield different values, since the 
map distortion is above the requested threshold given in ( [4T] ) 
and unphysical wiggles arise with the S-1 and S-3 schemes. At 
higher toroidal resolutions Aiptln <1/8 the N-3, S-1 and S-3 
schemes agree quite well, whereas the numerical perpendicu¬ 
lar diffusion for the N-1 scheme is large which causes a decay 
of the blob. At the highest toroidal resolution A^/ 27 r = 1/64 
the value drops slightly for the N-3 and S-1 schemes, which 
is an indication for numerical perpendicular diffusion. The 
convergence behaviour in dependence on poloidal resolution is 
shown in fig. [T^. It is apparent that the convergence for the 
N-1 scheme is slowest and could even not be achieved satis- 


Figure 14: Spatial resolution scan of diffusion of Gaussian blob 
= 0.88, Zg - -0.29, Wi? = wz = 8 • 10“^). Obtained maximum value on 
plane = tt at time t- 10.0 (For examples of snapshots see also fig. [12} in de¬ 
pendence on a) toroidal resolution for fixed poloidal resolution of h = 5 ■ 10“^ 
and b) in dependence on poloidal resolution for fixed toroidal resolution of 
Aip/27r= 1/32. 


3.2.3. Numerical diffusion 

Performing a resolution scan of numerical decay rates for 
zonal modes is again computationally elaborate. Moreover, the 
points per radial wavelength of a zonal mode varies within the 
poloidal planes due to flux expansion. Therefore, we consider 
only a single example and investigate if the numerical decay 
rates for realistic geometry correspond with the results from 
section 3.1.2 obtained in axial circular geometry. 

We consider a zonal mode r = 2, m = 0,n = 0 and choose 
a resolution of /z = 5 • 10 “"^ and Aipjln - 1 / 20 . The poloidal 
resolution of the mode varies between ~ 3 - 10 “^ at outboard 
midplane and kji 8 • 10“^ at the bottom. From the scaling 
( [47] ) and the numerical measurement given in fig. [TT^ we may 
expect decay exponents of: 


ynum ~ 


10 - 2 , 10-1 
10 - 4 , 10-2 
10-5,10-3 
10-^ 10-5 


for N-1, 
for N-3, 
for S-1, 
for S-3. 


(49) 
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We measured in the realistic geometry numerical decay rates of 
around: 


ynum 


8 • 10 -^ 
2 • 10 -^ 
4 • 10 -\ 
1 • 10 -^ 


which agrees with the estimates. 


for N-1, 
for N-3, 
for S-1, 
for S-3, 


(50) 


4. Conclusion and final remarks 

In the field line map approach field/fiux-aligned coordinates, 
which become singular on the separatrix/X-point, are avoided. 
The concept is based on a cylindrical grid, which is sparse in 
the toroidal direction, and a field line following discretisation 
for parallel operators to exploit the characteristic fiute mode 
property of the solutions. The discretisation of perpendicular 
operators is straight forward and simple, whereas in the dis¬ 
cretisation of parallel operators an interpolation (see equation 
or integration (see equation ( p^ ) is involved. Although it 
is not discussed here, there is no obvious reason, why the de¬ 
veloped methods could not easily be generalized also to three- 
dimensional equilibria, i.e. stellerators. 

Two discretisation schemes for the parallel diffusion opera¬ 
tor have been presented: A naive scheme, and a discretisation 
according to the support operator method, which conserves the 
self-adjointness property of the parallel diffusion operator on 
the discrete level. It has been shown in section 1231 with a two- 
dimensional model problem that the support schemes exhibit a 


lower (better scaling) numerical ’diffusion’ (see equations (34) 
and ([35])). The effects of a strongly distorted map on the sup¬ 
port scheme have been identified and verified. For interpolation 
based schemes (S-1 and S-3) a minimum toroidal resolution has 
to be supplied, to reduce the map distortion below the threshold 
given in expression gTJ. The problems with distorted maps can 
also be cured via the integration method for the parallel gradient 
(S-C). 

The numerical methods are implemented in the new code 
GRILLIX and extensive benchmarks, which show the validity 
of the field line map approach in general and GRILLIX in par¬ 
ticular, were presented in section Investigated were mainly 
the support and naive scheme each with a bilinear and a third 
order bipolynomial interpolation (N-1, N-3, S-1, S-3). Due to 
the interpolation, the convergence behaviour depends in general 
on the toroidal resolution A(f and also the poloidal resolution h. 


In agreement with the two-dimensional model (see section 2.5 
and I Appendix A. 3 ), the convergence rate with respect to the 
poloidal resolution is slowest for the N-1 scheme, whereas it 
is similar for the N-3, S-1 and S-3 scheme. Considering struc¬ 
tures with k\\ = 0, a general scaling for their numerical decay 
rate is suggested in equation ( [48] ). However, the presence of 
radial boundaries may introduce at high resolutions a break of 
this scaling for higher order interpolation methods. 

The scalings were derived only from a two-dimensional 
model and rigorous benchmarks have only been performed for 


axial circular geometry. However, the suggested results and 
scalings seem to be applicable also to realistic geometry, as 
long as the map distortion is below the proposed threshold (41), 
i.e. a sufficient toroidal resolution has to be supplied. As shown 
for example in section [333 numerical decay rates can be esti¬ 
mated also for general geometry from the scalings ( [48] ) and the 
numerical measurement shown in fig. [m 

From a practical point of view the S-3 scheme might be 
preferable, since it exhibits the lowest numerical diffusion 
among the presented methods. Even higher order interpolation 
methods (e.g. fifth order bipolynomial interpolation S-5) may 
be problematic, since their stencil can become large and they 
could cause unphysical oscillatory structures (Runge’s phe¬ 
nomenon) especially near the radial boundaries. Although the 
implementation of the S-C scheme might seem cumbersome, 
its advantage is that any spurious effects arising from map dis¬ 
tortion are absent. Future work might therefore address a com¬ 
bination of interpolation and integration, where the surface inte¬ 
gral in equation is executed by assuming interpolated poly¬ 
nomials as basis functions. 

In future work an application of the developed methods to a 
simple plasma turbulence model will be presented. 
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Appendix A. Additional model problems 

In this appendix additional two-dimensional model problems 
are discussed, where explicit expressions for the parallel diffu¬ 
sion operator are given. 

As in the model problem of section |2.5| v plays the role of 
a coordinate within poloidal planes and z the role of a peri¬ 
odic toroidal/axial coordinate. However, the magnetic field is 
assumed to be more general now: 

B = (A.l) 

with ^ Bo. Again the domain is spanned by a regular grid 
X/, Zk, with grid spacing h in the x-direction and grid spacing Az 
in the z-direction. From each grid point Xi,Zk the penetration 
points at the planes Zk±i are expressed via: 

xf = xi±h{nf+ft), (A.2) 
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where nf e Z, and 0 < /|- < 1 a dimensionless quantity which 
expresses the displacement of the penetration points with re¬ 
spect to the grid. For simplicity, we assume in the following 
that the lengths along field lines and the fiux box volumes are 
all equal, i.e. /S.Si = As and AVt = A'Vt = hAz. 


Appendix A.l. Inhomogeneous magnetic field (Interpolation) 

In section |2.5| the case for constant displacement has been 
discussed, i.e. fr=f = const.. We consider here the more 
general case where the displacement factors fi- vary. A linear 
interpolation is applied for the computation of the discrete par¬ 
allel gradient. 


The expression for the discrete parallel diffusion according to 
the naive scheme is given in fig. |A.15^ . The 'V and the 
discretisation for th e supp ort scheme yield in general different 


expressions. In fig. A. 15b the expression for is given 


and in fig. A. 15: for . The given expressions are valid 

only for nf = const. The final operator can be easily 

obtained as the average between both. 



Figure A. 15: Discrete parallel diffusion operator for non-tiomogeneous mag¬ 
netic field, a) Naive scheme b) support scheme discretisation 

^supp,^^ c)discretisation 0,7^’-. 


6 := aAzIBo 1 the penetration points are xf = xt (I ± e), and 
the displacement factors can be obtained as (see fig. |A.16 ): 

+ f-1, for: i < /, 

n^ =< 

1^0, for: i > /, 
l-e[i+(/-/)], for:i<I, 

+ (i - /)j, for; i> I ’ 

where xj = -hjl. Note that there is a jump present in nf at 
i = I. One may easily work out the expression for the dis¬ 
crete parallel diffusion operator. However, we give here only 
the result for a structure, which varies slowly in the v-direction, 
i.e. Ui^k = %• The result for the naive scheme is: 



(A.4) 

(A.5) 


(nnaive.~.\ _ + Uk-l + h+l 


(A.6) 


and for the support scheme: 


-f U]^-i -f Uk+i ^ ^ Uk-i 

As^ As 


hi) 


'2.Uk — 1 1 

A^2 


2A^ 

C Uk -\-1 U 

As 2As 


fori ^ 7,7+1, 

fori = 7,7+1. 
(A.7) 


The first term represents the discrete analogue of and the 
second term (V • b) V||W, since e/As ^ a/Bo V • b. Note that 
for / = 7,7+1 a small erroneous oscillation of size 6/4 arises due 
to the jump in nf. With the integration method for the discrete 
parallel gradient according to equation this oscillation does 
not arise (see section Appendix A.2 equation ( A.15| )). 

Such an oscillation is also observed with GRILLIX (S-1 and 
S-3). In practice, as this oscillation remains on the grid scale 
as long as \nf_^ ~ ^ it can also be cured by applying 

additionally a small perpendicular high order diffusion. 



An interesting case to consider is B^ — ax (This mimics the Figure A. 16: Displacement factors ff for model problem with diverging field 

situation at the X-point, where a convergent magnetic field in a lines, 
third dimension B^ = -ay would ensure V • B = 0). In the limit 
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Appendix A.2. Example for integration method 

We use the integration method (equation ( p^ ) for the dis¬ 
cretisation of the parallel gradient: 






(A.8) 


We consider the same example as described in [Appendix A.l 
The surface overlaps are (see fig. |A.17| ): 


aa;^, =h 


=h 


=h 


11 - 6 (/ - /) for i < /, 
1^1 - 6 (/ - / - 1) for / > /, 

jeQ - i 1) for i < /, 

1^0 for i > /, 

[o for/</, 

16 (/ - I) for i > /, 


Jl -6(/-/-r 1) for / </, 
AA • • =h < 

]l-e(i-I) 


AA-,_i =h- 




for i > I, 

0 for i < I, 

e(i -1 - 1) for / > I, 

e{i-1 - 1) for i < I, 

0 for i > 1, 


(A.9) 
(A.IO) 
(A.ll) 
(A.12) 
(A. 13) 
(A. 14) 





jl-e 


Figure A. 17: Surface overlaps AA^-^ for model problem with diverging field 
lines. 


One may again easily work out the expression for the discrete 
parallel diffusion operator, and we give again only the result for 
a structure which varies slowly in the v-direction i.e. 




-2% -f Uk-l + %+l _ € Uk+l - Uk-l ^ ^ 

a A.2 + A. 2A. 


In contrast to expression { A.l ) no oscillation around / = / is 
present. 


Appendix A. 3. Third order polynomial interpolation 

We consider again a homogeneous magnetic field, i.e. = 
ft - f - ^onst. but apply a third order polynomial inter¬ 
polation for the computation of the quantity at the penetration 
points. The problem, which contains at least one inner grid 
point is a 7 X 3 grid. The parallel gradient matrices are: 


= -(Q-f = 


/-I 


1 

A^ 


-1 


-1 


-1 


b a 
c b 
d c 
d 


-1 


b a 
c b 
d c 
d 


-1 


-1 


b a 
c b 
d c 
d 


-1 


-1 


b a 
c b 
d c 
d 


-1 


b a 


-1 


-1 


-1^ 


with: 


a = 


/7/-1) 

6 ’ 
(/-!)(/■ 


b = - 


/(/+!)(/- 2 ) 


2 ) 


/(/■ 


2 

1 )(/- 


2 ) 


(A. 16) 

z o 

The only inner point is / = 4, ^ = 2, which corresponds to 
the 11th row of the matrices Q-. The discrete parallel diffusion 
operator for the naive scheme is shown in fig. |A.18 f and for the 
support scheme in fig. [A.l8} ). Again, the stencil for the support 
scheme is bigger and if the penetration points coincide with grid 
points (/ = 0,1) both schemes yield again the standard second 
order finite difference expression. 

The action of the parallel gradient on a mode u = 
exp (fkxX + ik^z) yields: 


D" 


_ i.2 , J_i.4 A 2 _ 

ni ^ 12 II 


fe^)V(/-l)(/+!)(/-2) 

A^2 12 


+ 0 


(kxh, k^Azf 
A^2 


U, 


(A. 17) 


D 


supp 






(kxK kzAzf 

A^2 


u. 


(A.18) 


The error term k^^As^/l2 arises in both schemes and represents 
the discretisation error in the parallel direction. This error could 
be eliminated by using a fourth order finite difference expres¬ 
sion along magnetic field lines, where the stencil would cover 


also planes k±2. Although expression ( |A.18| ) might suggest that 
the error for the support scheme scales with respect to poloidal 
resolution like ik^h)^, it actually scales as the naive scheme also 
only like {kxh)^, since the parallel discretisation error is also de¬ 
pendent on {kxhf, i.e.: 


ktAs^ = 


(kx sin 6 kz cos 6)"^ {(fh)^ + Az^^ 


As^ 


(A. 19) 
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which contains also terms cx {kxhf /As^. Therefore, for k\\ ^ 0 
modes the schemes S-1, N-3 and S-3 exhibit the same scal¬ 
ing for the error with respect to poloidal resolution. However, 
only those terms where k\\ can not be factored out, represent a 
directional error and are responsible for numerical perpendicu¬ 
lar ’diffusion’. Therefore, for the naive scheme the numerical 
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support scheme like (kxh)^ /As^. 


a) 


As^ 


b) 


As^ 


+ 


+ 


w-w-m 






-r 

-r 

+ 


-r 


-r 




-r 


[-M+W-W\ 


'rM+W-m 






+ 


-r 

-r 

+ 


-r 

+ 


w-w-m 


+ 




[-fif+w-m 








|/V-lf^2)(f+iy36| 


+ 


-F 


-fif+w-m 


■^neticf.eia«n^ 


[8] H. R. Strauss, Edge localized mode simulations in divertor tokamaks, 
Physics of Plasmas 3 (1996) 4095. 

[9] B. D. Dudson, A. Allen, G. Breyiannis, E. Brugger, J. Buchanan, L. Easy, 
S. Farley, 1. Joseph, M. Kim, A. D. McGann, J. T. Omotani, M. V. Uman- 
sky, N. R. Walkden, T. Xia, X. Q. Xu, BOUT++: Recent and current 
developments. Journal of Plasma Physics 81 (2015) 365810104. 

[10] M. Ottaviani, An alternative approach to field-aligned coordinates for 
plasma turbulence simulations. Physics Letters A 375 (2011) 1677. 

[11] F. Hariri, M. Ottaviani, A flux-coordinate independent field-aligned ap¬ 
proach to plasma turbulence simulations. Computer Physics Communi¬ 
cations 184 (2013) 2419. 

[12] F. Hariri, P. Hill, M. Ottaviani, Y. Sarazin, The flux-coordinate inde¬ 
pendent approach applied to X-point geometries. Physics of Plasmas 21 
(2014) 082509. 

[13] S. Gunter, Q. Yu, J. Kruger, K. Lackner, Modelling of heat transport in 
magnetised plasmas using non-aligned coordinates. Journal of Computa¬ 
tional Physics 209 (2005) 354. 

[14] S. Gunter, K. Lackner, C. Tichmann, Finite element and higher order dif¬ 
ference formulations for modelling heat transport in magnetised plasmas. 
Journal of Computational Physics 226 (2007) 2306. 

[15] A. Stegmeir, D. Coster, O. Maj, K. Lackner, Numerical Methods for 3D 
Tokamak Simulations Using a Flux-Surface Independent Grid, Contribu¬ 
tions to Plasma Physics 54 (2014) 549. 

[16] P. J. M. Carthy, Analytic solutions to the Grad-Shafranov equation for 
tokamak equilibrium with dissimilar source functions. Physics of Plasmas 
6 (1999) 3554. 

[17] W. F. Ames, Numerical Methods for Partial Differential Equations, 3rd 
Edition, Academic Press, 1992. 

[18] M. Shashkov, S. Steinberg, Support-Operator Finite-Difference Algo¬ 
rithms for General Elliptic Problems, Journal of Computational Physics 
118 (1995) 131. 

[19] M. Shashkov, Conservative Finite-Difference Methods on General Grids, 
CRC Press, 1996. 

[20] D. Farina, R. Pozzoli, D. D. Ryutov, Effect of the magnetic field geometry 
on the flute-like perturbations near the X Point, Nuclear Fusion 33 (1993) 
1315. 

[21] T. T. Ribeiro, B. D. Scott, Conformal Tokamak Geometry for Turbulence 
Computations, IEEE Transactions on plasma science 38 (2010) 2159. 

[22] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical 
Recipes, 3rd Edition, Cambridge University Press, 2007. 

[23] J. M. Zerzan, Overlap: A Fortran Program for rapidly evaluating the Area 
of Overlap between two Polygons, Computers & Geosciences 15 (1989) 
1109. 


Figure A. 18: Discrete parallel diffusion operator for two-dimensional model 
problem with third order polynomial interpolation, a) Naive scheme, b) Support 
scheme. 


References 

[1] M. R. O’Brian, D. C. Robinson, Tokamak experiments, in: R. O. Dendy 
(Ed.), Plasma Physics: an Indroductory Course, Cambridge University 
Press, Cambridge, 1993. 

[2] P. C. Stangeby, G. M. McCracken, Plasma boundary phenomena in toka¬ 
maks, Nuclear Fusion 30 (1990) 1225. 

[3] K. Dietz, P.-H. Rebut, The ITER divertor, 15th lEEE/NPSS Symposium 
on Fusion Engineering 2 (1993) 604. 

[4] K. J. Dietz, S. Chiocchio, A. Antipenkov, G. Federici, G. Janeschitz, 
E. Martin, R. R. Parker, R. Tivey, Engineering and design aspects related 
to the development of the ITER divertor. Fusion Engineering and Design 
27 (1995) 96. 

[5] W. D. D’haeseleer, W. N. G. Hitchon, J. D. Callen, J. L. Sohet, Flux Co¬ 
ordinates and Magnetic Field Structure, Springer Series in Computational 
Physics, Springer, 1990. 

[6] E. Strumberger, S. Gunter, J. Hobirk, V. Igochine, D. Merkl, E. Schwarz, 
C. Tichmann, the ASDEX Upgrade Team, Numerical investigations of ax- 
isymmetric equilibria with current holes. Nuclear Fusion 44 (2004) 464. 

[7] N. Mattor, Coordinate system for use around an X-point, Physics of Plas¬ 
mas 2(1995) 594. 


15 


